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In gene mapping, it is common to test for association between the phenotype and the 
genotype at a large number of loci, i.e., the same response variable is used repeatedly 
to test a large number of non-independent and non-nested hypotheses. In many of these 
genetic problems, the underlying model is a mixed model consistent of one or very few 
major genes concurrently with a genetic background effect, usually thought as of polygenic 
nature and, consequently, modeled through a random effects term with a well-defined 
covariance structure dependent upon the kinship between individuals. Either because the 
interest lies only on the major genes or to simplify the analysis, it is habitual to drop the 
random effects term and use a simple linear regression model, sometimes complemented 
with testing via resampling as an attempt to minimize the consequences of this practice. 
Here, it is shown that dropping the random effects term has not only extreme negative 
effects on the control of the type I error rate, but it is also unlikely to be fixed by resampling 
because, whenever the mixed model is correct, this practice does not allow to meet some 
basic requirements of resampling in a gene mapping context. Furthermore, simulations 
show that the type I error rates when the random term is ignored can be unacceptably 
high. As an alternative, this paper introduces a new bootstrap procedure to handle the 
specific case of mapping by using recombinant congenic strains under a linear mixed 
model. A simulation study showed that the type I error rates of the proposed procedure 
are very close to the nominal ones, although they tend to be slightly inflated for larger 
values of the random effects variance. Overall, this paper illustrates the extent of the 
adverse consequences of ignoring random effects term due to polygenic factors while 
testing for genetic linkage and warns us of potential modeling issues whenever simple 
linear regression for a major gene yields multiple significant linkage peaks. 
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1. INTRODUCTION 

For more than four decades, linear mixed models have been 
used in a wide range of applications because of their conceptual 
simplicity and flexibility to accommodate correlated sources of 
variation as well as fixed regressors. A generic linear mixed model 
can be written as 

y = Xj8 + Zy + e (1) 

where X and Z are known incidence matrices, /? is a vector of 
unknown fixed regression coefficients, y is a vector of random 
effects, and e is the vector of errors. It is also common to assume 
that y and e are independent and both have null expectation and 
finite variances. In many situations, either intentionally or unin- 
tentionally, the statistical analysis is carried out ignoring the term 
Zy in the model. This practice, although recognized as inefficient, 
has been thought to be harmless whenever the interest resides 
solely on a subset of the regression coefficients with the remain- 
ing parameters of the model deemed as nuisance. This thought 
seems to be mostly based on the fact that fi° = (X'X) _1 X'y is 



still an unbiased and consistent estimator of j8. However, it is well 
known that ignoring Zy and using ordinary least squares, results 
in an estimator of Var(/J°) that is biased and inconsistent as well 
as non-independent of /J° (Dhymes, 1978). Of course, this will 
affect the distribution properties associated with f}° under nor- 
mality or, otherwise, the asymptotic properties of its distribution. 
It has been suggested that this problem can be mitigated if testing 
is done through resampling. However, the adverse consequences 
of dropping the random term from the mixed model is unlikely 
to be fixed by the use of resampling methods. In this paper, a spe- 
cific application to genetic mapping via recombinant congenic 
strains (RCS) of experimental animals is used to illustrate this. 
Briefly speaking, genetic mapping can be seen as a problem in 
which the association of one dependent variable (the phenotype) 
with a large number of potential explicative variables (the marker 
genotypes) is tested one-by-one or by taking a very small num- 
ber of markers at once. An RCS panel is a replicable mapping 
population for which animals within the same strain are con- 
sidered to be genetically identical and related to different degrees 
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with animals from other strains. Such an inter-strain relationship 
results in what is known as the genetic background effect and, 
whenever this effect is understood as the result of the addition of 
many components of minuscule effect, the inclusion of a random 
effects term in the model would be the natural way to account 
for it. 

A mouse panel of RCS is obtained by mating mice from two 
genetically distinct inbred strains (a donor strain and a recipient 
strain) followed by two or more rounds of backcrossing to the 
recipient strain and subsequent sister x brother mating without 
selection for particular markers or phenotypes for a minimum 
of 20 generations. The genetic resolution of the panel is con- 
trolled by the number of backcrossing rounds. Because of this 
construction, each strain of an RCS panel can be thought of as 
an inbred strain in which segments of random length from the 
genome of a recipient strain have been replaced with the corre- 
sponding segments from a donor strain. The main consequence 
of this breeding scheme is that non-linked genes controlling the 
same trait are separated and fixed in haplotypes of different 
strains, allowing the possibility of studying them individually. The 
standard RCS panel uses two backcross generations and, conse- 
quently, the total length of the segments from recipient strain 
constitute on average the 87.5% of the genome of each strain; 
the remaining 12.5% represents the total expected length of the 
replaced genome segments. Without loss of generality, this is the 
type of RCS considered in this paper. For a more comprehensive 
description of the RCS and their use in gene mapping see Demant 
and Hart (1986), Moen et al. (1992), and Fortin et al. (2001b, 
2007). Once the RCS panel have been established, the whole panel 
is genotyped to obtain full characterization of the genome of each 
strain. Each genotype data set can then be used for the anal- 
ysis of all individuals of the same strain; this is an important 
money-saving feature of the design since it does not require of re- 
genotyping each individual because, except for de novo mutations, 
all pups from the same strain are genetically identical. 

Although most mouse geneticists agree that RCS are a power- 
ful resource to map loci associated with complex traits, there is 
some disagreement on how to do the analysis. Originally, when 
the use of RCS for genetic mapping was proposed, the core idea 
was to look into the stain distribution pattern with respect to a 
phenotype of interest and identify the strain that exhibited the 
largest deviation from the other strains in the RCS panel and sub- 
sequently cross it with the recipient strain to obtain F\ and F2 
progenies to be analyzed by standard methods (Demant and Hart, 
1986; Fortin et al., 2001b). Two examples of the application of 
this approach are reported in Fortin et al. (2001a) and Mtillerova 
and Hozak (2004). The problem is that contrasting phenotypes 
from F\ mice versus the ones from the recipient strain will only 
be effective for dominant traits, while the power for additive traits 
will be diminished and lost completely for recessive traits. On the 
other hand, the analysis of the F% mice requires new genotyp- 
ing, which not only defeats the economic advantages of having 
developed RCS, but more importantly, because every F2 individ- 
ual has different genotype, this approach is not suited for complex 
quantitative traits when a single measurement may not be reli- 
able enough to determine the phenotype (Moen et al, 1992). 
Alternatively, there is a designs consisting of taking a sample of 



mice from each strain and analyzing the whole panel together. 
Although this approach does not require additional genotyping 
and has the potential for making more efficient use of the pheno- 
typic variation, also opens more room for analysis pitfalls if the 
proper model is not used. For example, Joober et al. (2002) uses 
a QTL mapping procedure equivalent to simple linear regression 
at the markers ignoring genetic background which, as pointed by 
Palmer and Airey (2003), it may result in false positive rates far 
in excess of the nominal value, even when Bonferroni corrections 
are used. Another common way to address the problem is to use 
strain averages as the phenotype and treat the panel of means as 
a backcross dataset for analysis purposes. This is essentially the 
"interval mapping" procedure proposed by Shao et al. (2010) and 
equivalent to the one used by Thifault et al. (2008). This approach 
may substantially reduce the power for RCS panels with reduced 
number of strains and it does not deal with the fact that the 
strains, related because their background, may not have the same 
kinship degree at genomic level and consequently the phenotype 
means may be not only non-independent but heteroscedastic, as 
well. Lee et al. (2006) and Camateros et al. (2010) extend the 
simple linear regression to account for the genetic background 
by adding a fixed factor ("background proportion" in the first 
paper; "background indicator" in the second). Although better 
than ignoring the background, from the genetics standpoint, it is 
difficult to justify the plausibility of a fixed effects model under 
the assumption that the background effect is the result of the 
additive action of many genes of minuscule effect. In fact, I argue 
that the natural way to model such a background effect consis- 
tent with the principles outlined by Fisher (1919) is through the 
inclusion of a random effects term in the model as implemented 
in Di Pietrantonio et al. (2010). In this paper, I describe in detail 
a procedure for the analysis of a quantitative trait locus (QTL) 
that models the genetic background (assumed to be of polygenic 
nature) as a random effect term and use this to show how the 
omission of such a term in the model leads to conclusions that 
are wrong and inconsistent with the data. 

2. MODELS 

2.1. THE NAIVE QTL MODEL FOR AN RCS PANEL 

In its simplest form, at each marker position m,m = 1, 2, . . . , M, 
the RCS/QTL model for the z'th individual, i = 1,2, ... ,n, can be 
written as 

y, 1 = M + lim Hm + e, (2) 

where y,- denotes the phenotype for the z'th individual, §,„ denotes 
the major locus effect associated with the mth marker, q\ m is the 
indicator of the BB genotype at the mth position which is deter- 
mined by the RCS data, and the e,s are a set of independent 
random variables with distribution A/"(0, cr 2 ) (AA and BB are 
the genotypes of the donor and recipient parental strain, respec- 
tively). Of course, under an oligogenic model, at most, a handful 
of £ m s should be different from zero. In fact, it is common prac- 
tice that at the first screening, the estimation is carried out by 
regression at each marker under the assumption of only one 
major gene. When the presumption of a dense enough genotyp- 
ing marker panel is not correct, procedures like modified interval 
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mapping can be used instead. Variations of the problem include 
conditioning on a given set of markers. The salient feature of this 
design is that, at the with marker position, one looks across the 
RCS panel and classifies each strain as either AA or BB, since 
under the model (Equation 2), this is the only source of genetic 
variation when estimating f m . However, this model ignores the 
fact that individuals from the same strain are genetically iden- 
tical (assuming no new mutation at the locus under scrutiny), 
and strains with the same ancestral background share large por- 
tions of their genome so that even without the involvement of a 
major gene, there is more likely to be reduced variation within 
strains. In a nutshell, regression mapping works by testing the 
association of the phenotype with the observed genotype at each 
marker location so that finding significant linkage at any position 
implies testing the M null hypotheses, f m = 0. Clearly, most of 
these hypotheses as well as their test statistics are not indepen- 
dent. This may lead to problems in the control of the type I error 
rate if multiple testing is not addressed properly. Another irregu- 
larity results from the fact that with a dense genotyping panel the 
number of tested hypotheses can by far exceed the sample size. 
Because of these considerations, p-value estimation by resampling 
of residuals has been seen as a plausible alternative. For this paper, 
the problem is addressed through bootstrap. 

2. 1. 1. Computation ofp-values 

The estimation of genome-wide corrected p-values by resam- 
pling requires that under the null hypothesis: (i) each resample 
is taken from an exchangeable distribution, (ii) the variation of 
the original sample is preserved through all resamples, and (iii) 
the genome-wide baseline for the test statistics at each position 
is the same. The first two requirements are standard for resam- 
pling in regression (Davison and Hinkley, 1997; Anderson and 
Ter Braak, 2003). The last requirement is imposed to ensure that 
the uncorrected p-values across the genome are comparable (this 
is particularly important when there are missing genotype data). 
One way to estimate corrected p-values is to select an ensemble 
of test statistics whose marginal distribution is the same when the 
model does not contain any major locus. 

Since under model (Equation 2) and the hypothesis of no 
major gene, the distribution ofy = (yi, yi y„)' is exchange- 
able, resampling from the raw observations will also preserve the 
variation through the pseudo-observations. This means that in 
the absence of non-genetic regressors or other non-oligogenic 
factors, resampling the raw phenotypes either by permutation 
or through bootstrap will produce similar results. Furthermore, 
under these premises, basic sampling and hypothesis testing prin- 
ciples indicate that a permutation based procedure will be more 
efficient and powerful. However, this is not necessarily the case 
when the premises are removed. Should the model also con- 
tain fixed non-genetic regressors, resampling from the leverage- 
adjusted residuals under the null hypothesis would be a procedure 
that approximates exchangeability while preserving the original 
variation of the data. However, under this situation, resam- 
pling from leverage-adjusted residuals results in a procedure 
with acceptable properties only in the bootstrap case (Davison 
and Hinkley, 1997), while this is not longer guaranteed when 
resampling via permutation. The main issue is that sampling 



without replacement magnifies the effects of modest departures 
from exchangeability. Then, permuting leverage-adjusted residu- 
als may not be good enough (even worst, it may not be valid) 
and we would require of a much more elaborate and computer 
intensive procedure to obtain residuals guaranteed to be at least 
weakly exchangeable so that permutation works properly (see, 
for example, Kherad-Pajouh and Renaud, 2010). To complete 
the requirements listed above regarding the possibility of miss- 
ing genotypes, we propose to use the test statistic defined by the 
expression 

z m = t m 1 - - — 1 + - — where t m = - — (3) 

\ 4v m J \ 2v m J ff| m 

and £ m is the ordinary least squares estimate of £ m , m = 
1, 2, . . . , M, i.e., z m is just t m , our familiar f-statistic with v m 
degrees of freedom, transformed into a z-score (v m may vary 
slightly from marker to marker due to missing data). Another 
option would be a modified f-statistic t' m in which the mth esti- 
mate of variance s 2 m used to compute a| is replaced by sjj, the 

estimate under the null hypothesis. With no missing genotypes 
the use of any of z m , t' m , and t m would yield approximately the 
same p-value estimates. 

2. 1.2. Bootstrap procedure for simple linear regression at the 
markers 

The following bootstrap procedure computes the genome-wide 
corrected p-values for model (Equation 2) with the test statistic 
(Equation 3): 

STEP 1. At each marker position, m, fit the simple linear regres- 
sion at the markers model (Equation 2), use (Equation 3) 
to compute the test statistic z,„, and obtain the genome- 
wide set of statistics Zm = {z m , m= 1,2,..., M}. Also, 
set the genome-wide acceptance count vector to zero. 

Step 2. Sample with replacement from the raw vector of phe- 
notypes, y € W, to obtain y* e W, a bootstrapped full 
replica of y, and use this vector to compute z^^ = 
max {z* m }, where z*,, m = 1, 2, . . . , M, is the test statis- 
tic at the mth locus, computed by using y*, the vector of 
the pseudo-observations, instead of the original vector of 
phenotypes. 

Step 3. For each z m in Zm, if z m < z^ ax , add a unit to the mth 
entry of the acceptance count vector. 

Step 4. Repeat steps 1 and 2 R times and then compute the esti- 
mate of the vector of p-values by dividing the acceptance 
count vector by _R. 

This resampling scheme can be seen as an adaptation of a reg- 
ular regression residuals bootstrapping procedure (Davison and 
Hinkley, 1997), coupled with Roy's union-intersection principle 
(Roy, 1953) to control for the genome-wide type I error rate. 
When applied to the analysis of the RCS panel, this procedure is 
valid when there is only one observation per strain or when the 
within-strain variation is negligible. Otherwise, a random term 
in the model has been neglected and, regardless of | m being an 
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unbiased estimator of f m , the exchangeability requirement can- 
not be met and the most likely consequence would be an inflated 
type I error rate. In fact, as per arguments given by Churchill 
and Doerge (1994) and Churchill and Doerge (2008), this state- 
ment is correct not only for the bootstrap and RCS, but also for 
permutation test procedures applied to any study design involv- 
ing replicable mapping populations because, as for bootstrap, the 
Fisher (1935) principle of permutation also relies on exchange- 
ability. For simple experimental designs such as an intercross or a 
backcross mating, the individual units can safely be assumed to be 
exchangeable. However, it would be wrong to assume exchange- 
ability for more complicated designs, like advanced intercross, 
heterogeneous stocks and RCS. 

2.2. THE QTL MIXED MODEL FOR AN RCS PANEL 

The previous simple linear model (Equation 2) generalizes to a 
model of the form: 



y = XP + Zy+ q m £ m + e 



(4) 



where y represents the phenotype vector, q m is a vector with 
each entry being an indicator variable of the genotype BB at the 
marker position m with £ m being its associated effect (major gene 
effect), y is a random effects vector associated with the genetic 
background with E(y) = 0 and Var(y) = a 2 Ai, with a 2 > 0 
and Ai, a positive-definite matrix, both assumed to be constant, 
although unknown, X is a matrix of fixed covariates and its corre- 
sponding parameter vector /?, e is a vector of independent and 
identically distributed random variables representing the error 
term with E(e) =0 and Var(e) = a 2 1. Up to a multiplicative con- 
stant, Ai is a function of the length of the segments identical 
by descent shared amongst strains. For an established RCS panel 
there are only two possible identity states between pairs of strains 
at a given locus: either (i) all four alleles are identical by descent 
( Ai is the matrix holding the pairwise probabilities for this state), 
or (ii) the strains have different allelic forms and thus identical by 
descent only amongst themselves. So an estimator of Ai with "a 
high degree of precision" can be reached. Such an estimator uses 
only genomic information and does not involve y, so when esti- 
mating the parameters, one can assume that Ai is given. Another 
option is to take the entries of A i as the expected value of the 
proportion of the genome shared identical by descent between 
the respective strains under the RCS panel construction described 
above, i.e., 



1 i£i = j 

if i and j have the same background (5) 



15 
16 

if i and j have different backgrounds. 



This option, although not the most efficient, does capture the 
main features of the design and yields a variance structure for the 
random effects vector that can be exploited in the implementa- 
tion of the resampling algorithm. For example, if all the strains 
in the panel under scrutiny have the same background and the 
simplified expectation-based A i is used, then the distribution of 
the vector of random effects is exchangeable. Nonetheless, replac- 
ing a genomic-based Ai estimate by its theoretical expectation 



(Equation 5) implies ignoring important information regarding 
the correlation of the additive polygenic effects associated to the 
genetic background. 

2.2.1. Estimation 

The estimation for the mixed linear model has been extensively 
discussed in the literature (Harville, 1977; Henderson, 1986). 
Here we develop an application of these standard methods to the 
RCS design. Without loss of generality, let us consider the linear 
mixed model (Equation 1) with Var(y) 
o 2 l. Thus 



az Ai and Var(e) 



E(y) = X/? and Var(y) = a 2 (ZGZ' + i) = a 2 £ 

where G = k&\ and A. = i.e., k represents the signal-to-noise 
ratio. Under the assumption of no major gene and only polygenic 
background, k is related to the heritability coefficient. When G is 
known, the best linear unbiased estimator of /? and the best linear 
unbiased predictor of y (also known as a shrinkage estimator) can 
be written as 

p = (W'WrWv and y = GZ'Z z (v - Wj8), 

respectively, where W = T, ~ z X and v = E ~ z y. Also 
1 



a 2 = 



N - rank(W) 



(v-WjS)'(v- W/J) 



g.2 _ _ 

7 ~ rank(G) 



{y'G- l y +o- 2 tr(G" 1 C)) 



with 



C = (Z'MZ + G -1 ) -1 and M = I — X(X'X) X'. 

Notice that the previous expressions cannot be computed unless 
the signal-to-noise ratio, k, is known. A situation of a more prac- 
tical interest is an iterative procedure on which k is replaced by 
its estimate and, once that the estimates of a 2 and a 2 have been 
updated, a refinement of the estimate of k is obtained and so on. 
This iterative procedure will result in a fi and y that are no longer 
linear, nonetheless, they preserve most of the desirable properties 
present in their linear counterpart (Jiang, 1998). 

2.2.2. Mixed model resampling scheme 

Let us now focus our attention toward a resampling scheme 
appropriate for RCS data under a mixed model. By now, it is 
obvious that the bootstrap procedure described in the previous 
section will not work for the mixed model (Equation 4). A crude 
extension to this procedure would consist of computing 

e = y - X/? - Zy 

and resampling from y and e to obtain y* and e* so that the 
pseudo-observation y* could be recovered as 

y* =X/j + Zy* + e*. 
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However, it is straightforward to see that these residuals are not 
exchangeable and they are biased toward zero. Thus, they may 
not adequately represent the hypothesis tested nor reflect the true 
variation of the model. 

Alternatively, note that when p and k are known, it follows 
from the model under the null hypothesis that E(v) = W/J and 
Var(v) = a 2 I which implies that the distribution of the vec- 
tor of residuals, e = v — Wj8, is exchangeable. This suggests the 
following residuals resampling scheme: 

(1) given X and (S obtained under the mixed model without a 
major gene, i.e., under the null hypothesis, compute E, W by 
replacing k with k and A i with its genomic-based estimate; 
then, obtain the leverage-adjusted residuals 

€ = D(S"Jy-Wj8) 

where D is a diagonal matrix with each of the non-zero 
elements given by (1 — and ha is the rth leverage 

coefficient; 

(2) with replacement, resample from e e R" to obtain e* G W, 
its bootstrapped replica, and construct the vector of pseudo- 
observations as 

v* = Wj8 + e*. 

If instead of a bootstrap procedure based on leverage-adjusted 
residuals we want to use a residuals-based permutation proce- 
dure, then we need to extend the method of Kherad-Pajouh and 
Renaud (2010) to get weak exchangeability of residuals. However, 
when k is estimated from the data, such an extension is not possi- 
ble and we would have to rely on approximations. More research 
is needed to explore this direction. 

Outside of a genetics context, there is a number of permuta- 
tion and bootstrap procedures for mixed models whose objective 
is testing the components of variance (for example, Fitzmaurice 
et al., 2007; Sinha, 2009; Lee and Braun, 2012; Samuh et al., 
2012). However, they cannot be applied in our case because we 
are interested in the regression coefficients (or a subset of them) 
and the variance of the random effects is just nuisance parameter. 
Incidentally, when testing the components of variance, bootstrap 
has the edge over most permutation procedures (Samuh et al., 
2012). 

2.2.3. Bootstrap procedure for the mixed linear model 

According to the foregoing argument, generalization to the previ- 
ous bootstrap procedure to compute the genome-wide corrected 
p-values for the mixed model (Equation 4) goes as follows: 

Step 0. Compute Ai from the genotype data of the RCS panel, 
and under the null hypothesis, obtain k, p, E, W and e 
as described in (i) above. 

Step 1. At each marker position, m, fit the model 

*-(* r S(£) +e ' (6) 

Of course, this model is equivalent to model 
(Equation 4), the RCS/QTL mixed model, with k 



replaced by k. Compute the model parameter estimates 
with the outlined mixed model procedure as well as the 
test statistic set Z = {z m , m = 1,2, ... , M] by using 
Equations (6) and (3); set the acceptance count vector to 
zero. 

Step 2. Draw a pseudo-observation v* by using the proposed 
resampling scheme in (ii) above and fit the major gene 
model in model (Equation 6) with v replaced by v* to 
obtain the set of bootstrapped test statistics {z* m } and its 
associated critical value z* „ = max \z* }. 

lilil.V 1 if I* 

Step 3. For each z m in Z, if z m < zj^ax' add a unit to the mth 

entry of the acceptance count vector. 
Step 4. Repeat steps 2 and 3 _R times and compute the p-value 

estimates by dividing the acceptance count vector by -R. 

To my knowledge, this bootstrap procedure for the analyzing a 
panel of RCS has not been proposed before Di Pietrantonio et al. 
(2010) and this paper contains the first detailed derivation and 
study of its properties. In fact, the resampling methods (mostly 
conditional permutation) applied to analyze RCS have not used 
mixed models, but consider the strain effect as fixed which is 
inconsistent with the hypothesis of a genetic background of poly- 
genic nature or discard information by using only the estimated 
strain means (for example, Gill and Boyle, 2005; Thifault et al., 
2008; Camateros et al, 2010). 

3. RESULTS 

One straightforward way to show the effect of ignoring the ran- 
dom effects term in a mixed model is by simulation. The idea 
is to generate a dataset from a model that includes a random 
term for genetic background and noise, but is free of any major 
locus. Then compare thep-value profiles (actually, — log 10 p pro- 
files) obtained by the use of the naive model (Equation 2) as well 
as the mixed model (Equation 4). For this simulation study, the 
genotypes of an RCS panel of 36 strains that were described in 
Fortin et al. (2001b) were used. The panel originally had 37 lines 
and 625 microsatellite markers; since then, one line has died out 
and six markers were removed for reliability reasons. Although a 
much larger set of single nucleotide polymorphism markers for 
this RCS panel is also available, I think that this set of 619 mark- 
ers is enough to show the harmful effects of fitting the wrong 
model on the inference. Of course, more markers will only exac- 
erbate the problem. For this simulation experiment, six different 
values for the signal-to-noise ratio parameter k were chosen (0, 

4, |, j, 1, and 2). Under a standard additive polygenic model, 
i.e., a model without major genes, the signal-to-noise parameter 
is a function of the heritability coefficient (the chosen values cor- 
respond to the heritability proportions of 0, i, ^, i, j, and I, 
respectively). In every simulation run, a sample of seven individ- 
uals from each strain was simulated under the assumption of no 
major gene, i.e., under model (Equation 4) with £,„ = 0 for all 
markers, m = 1,2, ... ,M. The value of a 2 was fixed for all sim- 
ulations to 1.175, while X/J was fixed as a vector with 7 in all its 
entries. Simulations for each value of k were run 1000 times and 
both methodologies, the mixed model as well as the bootstrapped 
naive regression at the markers were applied to the simulated 
datasets with 10, 000 as the number of resamples for every dataset. 
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In gene mapping studies, a significant peak is defined as the most 
extreme point of a region beyond the p-value threshold according 
to some pre-specified genome-wide type I error rate (Churchill 
and Doerge, 1994). For this study, we use a value of 0.01 or equiv- 
alently, a threshold value of 2 on a — log 10 p-scale. Tables 1-3 
summarize the results of these simulations. As expected, when- 
ever there is not a polygenic term in the model (i.e., A = 0), both 
methodologies produce identical results. However, the picture 
changes when A > 0. In this case, it is quite obvious that ignor- 
ing the random effects term has pernicious consequences even for 
modest levels of A, the signal-to-noise ratio, while the proposed 
mixed model method keeps the genome-wide type I error rate 
relatively close to the nominal value. However, the empirical type 
I error rates obtained by the proposed procedure seem to increase 
slightly with A (Table 3). This phenomenon maybe due to the fact 
that the makers used for mapping purposes are also used to esti- 
mate the probability of identity by descent between strains and, to 
a lesser extent, the fact that the the bootstrap procedure is based 
on residuals computed with X and /? estimated from the same 
data. Nonetheless, the moral of this exercise is that whenever sim- 
ple regression of a major gene model produces many significant 
peaks, a warning flag about the model validity should be raised. 



Table 1 | Percentage of declared significant peaks with a bootstrap 
genome-wide adjusted significance level of 0.01 when the proposed 
mixed model methodology is used. 
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Table 2 | Percentage of declared significant peaks with a bootstrap 
genome-wide adjusted significance level of 0.01 when a naive 
regression at the markers is used. 
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The histogram of a typical dataset obtained by simulation from 
a model with polygenic effects only would look like the one shown 
in Figure 1. Nonetheless, for this histogram I chose a dataset 
for which simple linear regression produces a very large num- 
ber of significant peaks. If a major locus were at play, one would 
expect to have a well-defined bimodal distribution, so this his- 
togram seems consistent with the generating model of no major 
gene. However, when we look into the p-value profiles obtained 
through the model that ignores the genetic background term, 
instead of profiles consistent with the model we will have some- 
thing extreme as shown by dashed lines in Figure 2. According to 
the profiles on this figure, one might conclude that all chromo- 
somes have at least one significant peak, fact that does not appear 
to be supported by the histogram of the data, and more conclu- 
sively, this is in conflict with the generating model. If anything, it 
can be argued that the data distribution may seem a bit skewed, 
but one may expect that estimation of p-values via bootstrapping 
of residuals should not be too sensitive to this. Of course, as for 
bi-modality, skewness may also be caused by a mixture of distri- 
butions. However, a very strong peak, as any of the ones spotted 
on every chromosome, is difficult to conceive without a conspic- 
uous bimodal distribution. Even with the use of robust regression 
estimates instead of the obtained by regular least squares to min- 
imize the potential impact of outliers on the estimation, these 
profiles change very little (data not shown). When the missing 
random effects term is introduced into the model (solid blue 



Table 3 | Empirical genome-wide type 1 error rates obtained via 
bootstrap in the simulation study (0.01 is the nominal value and the 
number of simulated datasets for each X is 1000). 
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FIGURE 1 | Typical histogram of simulated data. The p-value profiles of 
the data on this histogram were computed and plotted in Figure 2. 
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FIGURE 2 | Bootstrap genome-wide corrected p-value profiles. Dashed line for naive model (Equation 2) and solid line (sometimes hardly distinguishable 
from the x-axis line) for the mixed model (Equation 6). Note that both profiles have been corrected for multiple testing. 



lines in Figure 2), p-value profiles become consistent with the 
generating model. Repetition of this exercise on any other simu- 
lated datasets yields similar results, although the specific resulting 
profiles most likely are not be the same. 

4. DISCUSSION 

This paper proposes a bootstrapping procedure to estimate thep- 
values under a mixed model applied to gene mapping when RCS 
are used. The method can be easily adapted for other replicable 
mapping population/designs. This procedure is a generalization 
of the linear regression bootstrap of residuals coupled with the 
union-intersection principle aimed to control the genome-wide 
type I error rate. A simulation study with different values of the 
signal-to-noise ratio unequivocally shows that when a panel of 



RCS is used for mapping, ignoring one random effects term in a 
mixed linear model can have pernicious consequences, resulting 
in inflated type I error rates and leading to the declaration of sig- 
nificant linkage peaks were no such peaks should be found. The 
simulation study also shows that the proposed bootstrap proce- 
dure seems to produce slightly inflated type I error rates as the 
signal-to-noise ratio increases. This problem is likely due to the 
fact that the markers used for mapping are also used to estimate 
the length of the segments shared identical by descent but also 
it can be associated with a stronger departure from exchange- 
ability as the ratio increases. In any case, the problem deserves 
further scrutiny. The proposed bootstrap procedure for mixed 
models is quite general and can easily be adapted to non-genetic 
problems. 
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